      SUBROUTINE PARNSET
C------------------------------------------------------------------------
C
C
C------------------------------------------------------------------------
      INCLUDE 'AVL.INC'
C
C---- set parameter unit name
      DO IP = 1, IPTOT
        PARUNCH(IP) = ' '
      ENDDO
C
      PARUNCH(IPALFA) = 'deg'
      PARUNCH(IPBETA) = 'deg'
      PARUNCH(IPPHI)  = 'deg'
      PARUNCH(IPTHE)  = 'deg'
      PARUNCH(IPPSI)  = 'deg'
      PARUNCH(IPVEE)  = UNCHV
      PARUNCH(IPRHO)  = UNCHD
      PARUNCH(IPGEE)  = UNCHA
      PARUNCH(IPRAD)  = UNCHL
      PARUNCH(IPXCG)  = UNCHL
      PARUNCH(IPYCG)  = UNCHL
      PARUNCH(IPZCG)  = UNCHL
      PARUNCH(IPMASS) = UNCHM
      PARUNCH(IPIXX)  = UNCHI
      PARUNCH(IPIYY)  = UNCHI
      PARUNCH(IPIZZ)  = UNCHI
      PARUNCH(IPIXY)  = UNCHI
      PARUNCH(IPIYZ)  = UNCHI
      PARUNCH(IPIZX)  = UNCHI
C
      RETURN
      END ! PARNSET



      SUBROUTINE DEFINI
C------------------------------------------------------------------------
C
C
C------------------------------------------------------------------------
      INCLUDE 'AVL.INC'
C
C---- Does this really need to be done at runtime?
      PI = 4.0*ATAN(1.0)
      DTR = PI/180.0
C
C---- logical units
      LUINP = 4   ! configuration file
      LURUN = 7   ! run case file
      LUMAS = 8   ! mass file
      LUPRM = 9   ! parameter file
      LUOUT = 19  ! output dump file
      LUSTD = 20  ! stability derivative dump file
      LUSYS = 22  ! dynamic system matrix dump file
C
C---- flag for forces in standard NASA stability axes (as in Etkin)
      LNASA_SA  = .TRUE.
C
C---- flag for rotations defined in stability axes or body axes
      LSA_RATES = .TRUE.
C
      LPTOT   = .TRUE.
      LPSURF  = .FALSE.
      LPSTRP  = .FALSE.
      LPELE   = .FALSE.
      LPHINGE = .FALSE.
      LPDERIV = .FALSE.
C
      LGEO  = .FALSE.
      LENC  = .FALSE.
C
      LAIC  = .FALSE.
      LSRD  = .FALSE.
      LVEL  = .FALSE.
      LSOL  = .FALSE.
      LSEN  = .FALSE.
C
      LVISC    = .TRUE.
      LBFORCE  = .TRUE.
      LTRFORCE = .TRUE.
C
      LMWAIT = .FALSE.
C
      MATSYM = 0
      NITMAX = 20
C
      SAXFR = 0.25  ! x/c location of spanwise axis for Vperp definition
C
      VRCORE = 0.25   ! vortex core radius / vortex span
      SRCORE = 0.75   ! source core radius / body radius
C
C---- dafault basic units
      UNITL = 1.
      UNITM = 1.
      UNITT = 1.
      UNCHL = 'Lunit'
      UNCHM = 'Munit'
      UNCHT = 'Tunit'
      NUL = 5
      NUM = 5
      NUT = 5
C
C---- set corresponding derived units
      CALL UNITSET
C
C---- default air density and grav. accel.
      RHO0 = 1.0
      GEE0 = 1.0
C
C---- no eigenvalue reference data yet
      FEVDEF = ' '
      DO IR = 1, NRMAX
        NEIGENDAT(IR) = 0
      ENDDO
C
C---- no run cases defined yet
      NRUN = 0
      IRUN = 1
C
C---- number of valid time levels stored
      NTLEV = 0
C
C---- default time step, and number of time steps to take
      DELTAT = 0.0
      NTSTEPS = 0
C
      RETURN
      END ! DEFINI



      SUBROUTINE PARSET
C------------------------------------------------------------------------
C
C
C------------------------------------------------------------------------
      INCLUDE 'AVL.INC'
C
C---- variable names
      VARNAM(IVALFA) = 'alpha '
      VARNAM(IVBETA) = 'beta  '
      VARNAM(IVROTX) = 'pb/2V '
      VARNAM(IVROTY) = 'qc/2V '
      VARNAM(IVROTZ) = 'rb/2V '
C
C---- variable selection keys
      VARKEY(IVALFA) = 'A lpha'
      VARKEY(IVBETA) = 'B eta'
      VARKEY(IVROTX) = 'R oll  rate'
      VARKEY(IVROTY) = 'P itch rate'
      VARKEY(IVROTZ) = 'Y aw   rate'
C
C---- constraint names
CCC                     123456789012
      CONNAM(ICALFA) = 'alpha '
      CONNAM(ICBETA) = 'beta  '
      CONNAM(ICROTX) = 'pb/2V '
      CONNAM(ICROTY) = 'qc/2V '
      CONNAM(ICROTZ) = 'rb/2V '
      CONNAM(ICCL  ) = 'CL    '
      CONNAM(ICCY  ) = 'CY    '
      CONNAM(ICMOMX) = 'Cl roll mom'
      CONNAM(ICMOMY) = 'Cm pitchmom'
      CONNAM(ICMOMZ) = 'Cn yaw  mom'
C
C---- constraint selection keys
      CONKEY(ICALFA) = 'A '
      CONKEY(ICBETA) = 'B '
      CONKEY(ICROTX) = 'R '
      CONKEY(ICROTY) = 'P '
      CONKEY(ICROTZ) = 'Y '
      CONKEY(ICCL  ) = 'C '
      CONKEY(ICCY  ) = 'S '
      CONKEY(ICMOMX) = 'RM'
      CONKEY(ICMOMY) = 'PM'
      CONKEY(ICMOMZ) = 'YM'
C
C------------------------------------------------------------------------
      IZERO = ICHAR('0')
C
C---- add control variables, direct constraints
      DO N = 1, NCONTROL
        ITEN = N/10
        IONE = N - 10*ITEN
C
C------ assign slots in variable ond constraint lists
        IV = IVTOT + N
        IC = ICTOT + N
        VARNAM(IV) = DNAME(N)
        CONNAM(IC) = DNAME(N)
        IF(ITEN.EQ.0) THEN
         VARKEY(IV) = 'D' // CHAR(IZERO+IONE) // ' '
     &             // ' ' // DNAME(N)(1:8)
         CONKEY(IC) = 'D' // CHAR(IZERO+IONE)
        ELSE
         VARKEY(IV) = 'D' // CHAR(IZERO+ITEN) // CHAR(IZERO+IONE)
     &             // ' ' // DNAME(N)(1:8)
         CONKEY(IC) = 'D' // CHAR(IZERO+ITEN) // CHAR(IZERO+IONE)
        ENDIF
C
        LCONDEF(N) = .TRUE.
      ENDDO
C
C---- default design-variable flags, names
      DO N = 1, NDESIGN
        LDESDEF(N) = .TRUE.
      ENDDO
C
C---- total number of variables, constraints
      NVTOT = IVTOT + NCONTROL
      NCTOT = ICTOT + NCONTROL
C
C---- run-case parameter names
      PARNAM(IPALFA) = 'alpha    '
      PARNAM(IPBETA) = 'beta     '
      PARNAM(IPROTX) = 'pb/2V    '
      PARNAM(IPROTY) = 'qc/2V    '
      PARNAM(IPROTZ) = 'rb/2V    '
      PARNAM(IPCL )  = 'CL       '
      PARNAM(IPCD0)  = 'CDo      '
      PARNAM(IPPHI)  = 'bank     '
      PARNAM(IPTHE)  = 'elevation'
      PARNAM(IPPSI)  = 'heading  '
      PARNAM(IPMACH) = 'Mach     '
      PARNAM(IPVEE)  = 'velocity '
      PARNAM(IPRHO)  = 'density  '
      PARNAM(IPGEE)  = 'grav.acc.'
      PARNAM(IPRAD)  = 'turn_rad.'
      PARNAM(IPFAC)  = 'load_fac.'
      PARNAM(IPXCG)  = 'X_cg     '
      PARNAM(IPYCG)  = 'Y_cg     '
      PARNAM(IPZCG)  = 'Z_cg     '
      PARNAM(IPMASS) = 'mass     '
      PARNAM(IPIXX)  = 'Ixx      '
      PARNAM(IPIYY)  = 'Iyy      '
      PARNAM(IPIZZ)  = 'Izz      '
      PARNAM(IPIXY)  = 'Ixy      '
      PARNAM(IPIYZ)  = 'Iyz      '
      PARNAM(IPIZX)  = 'Izx      '
      PARNAM(IPCLA)  = 'visc CL_a'
      PARNAM(IPCLU)  = 'visc CL_u'
      PARNAM(IPCMA)  = 'visc CM_a'
      PARNAM(IPCMU)  = 'visc CM_u'
C
C---- total number of parameters
      NPTOT = IPTOT
C
C---- set default parameter unit names
      CALL PARNSET
C
      RETURN
      END ! PARSET





      SUBROUTINE VARINI
C------------------------------------------------------------------------
C
C
C------------------------------------------------------------------------
      INCLUDE 'AVL.INC'
C
C---- initialize state
      ALFA = 0.
      BETA = 0.
      WROT(1) = 0.
      WROT(2) = 0.
      WROT(3) = 0.
C
      DO N = 1, NCONTROL
        DELCON(N) = 0.0
      ENDDO
C
      DO N = 1, NDESIGN
        DELDES(N) = 0.0
      ENDDO
      LSOL = .FALSE.
C
      RETURN
      END ! VARINI






      SUBROUTINE RUNINI
C------------------------------------------------------------------------
C
C
C------------------------------------------------------------------------
      INCLUDE 'AVL.INC'
C
      WRITE(*,*)
      WRITE(*,*) 'Initializing run cases...'
C
C---- go over all run cases
      DO IR = 1, NRMAX
C------ index of default constraint for each variable
        ICON(IVALFA,IR) = ICALFA
        ICON(IVBETA,IR) = ICBETA
        ICON(IVROTX,IR) = ICROTX
        ICON(IVROTY,IR) = ICROTY
        ICON(IVROTZ,IR) = ICROTZ
C
C------ default constraint values
        DO IC = 1, ICTOT
          CONVAL(IC,IR) = 0.
        ENDDO
C
C------ default run case titles
        RTITLE(IR) = ' -unnamed- '
C
C------ default dimensional run case parameters
        DO IP = 1, NPTOT
          PARVAL(IP,IR) = 0.
        ENDDO
        PARVAL(IPGEE,IR) = GEE0
        PARVAL(IPRHO,IR) = RHO0
C
C------ default CG location is the input reference location
        PARVAL(IPXCG,IR) = XYZREF0(1)
        PARVAL(IPYCG,IR) = XYZREF0(2)
        PARVAL(IPZCG,IR) = XYZREF0(3)
C
        PARVAL(IPMASS,IR) = RMASS0
        PARVAL(IPIXX,IR) = RINER0(1,1)
        PARVAL(IPIYY,IR) = RINER0(2,2)
        PARVAL(IPIZZ,IR) = RINER0(3,3)
        PARVAL(IPIXY,IR) = RINER0(1,2)
        PARVAL(IPIYZ,IR) = RINER0(2,3)
        PARVAL(IPIZX,IR) = RINER0(3,1)
C
        PARVAL(IPCD0,IR) = CDREF0
C
        PARVAL(IPCLA,IR) = DCL_A0
        PARVAL(IPCLU,IR) = DCL_U0
        PARVAL(IPCMA,IR) = DCM_A0
        PARVAL(IPCMU,IR) = DCM_U0
C
        ITRIM(IR) = 0
        NEIGEN(IR) = 0
      ENDDO
C
C---- add control variables, direct constraints
      DO N = 1, NDMAX
        IV = IVTOT + N
        IC = ICTOT + N
        DO IR = 1, NRMAX
          ICON(IV,IR) = IC
          CONVAL(IC,IR) = 0.
        ENDDO
      ENDDO
C
C---- default number of run cases
      IRUN = 1
      NRUN = 1
C
C---- all run cases are targets for eigenmode calculation
      IRUNE = 0
C
C---- first run case is default for time march initial state
      IRUNT = 1
C
      RETURN
      END ! RUNINI



      SUBROUTINE RUNGET(LU,FNAME,ERROR)
C-------------------------------------------------
C     Reads run case file into run case arrays
C-------------------------------------------------
      INCLUDE 'AVL.INC'
      CHARACTER*(*) FNAME
      LOGICAL ERROR
C
      CHARACTER*80 LINE
      CHARACTER*12 VARN, CONN
      CHARACTER*8  PARN
C
      OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=90)
      ILINE = 0
C
      IR = 0
C
C==============================================================
C---- start line-reading loop
 10   CONTINUE
C
      READ(LU,1000,END=50) LINE
 1000 FORMAT(A)
      ILINE = ILINE + 1
C
      KCOL = INDEX(LINE,':' )
      KARR = INDEX(LINE,'->')
      KEQU = INDEX(LINE,'=' )
      IF(KCOL.NE.0) THEN
C----- start of new run case
       READ(LINE(KCOL-3:KCOL-1),*,ERR=80) IR
C
       IF(IR.LT.1 .OR. IR.GT.NRMAX) THEN
        WRITE(*,*) 'RUNGET:  Run case array limit NRMAX exceeded:', IR
        IR = 0
        GO TO 10
       ENDIF
C
       NRUN = MAX(NRUN,IR)
C
       RTITLE(IR) = LINE(KCOL+1:80)
       CALL STRIP(RTITLE(IR),NRT)
C
      ELSEIF(IR.EQ.0) THEN
C----- keep ignoring lines if valid run case index is not set
       GO TO 10
C
      ELSEIF(KARR.NE.0 .AND. KEQU.NE.0) THEN
C----- variable/constraint declaration line
       VARN = LINE(1:KARR-1)
       CONN = LINE(KARR+2:KEQU-1)
       CALL STRIP(VARN,NVARN)
       CALL STRIP(CONN,NCONN)
C
       DO IV = 1, NVTOT
         IF(INDEX(VARNAM(IV),VARN(1:NVARN)).NE.0) GO TO 20
       ENDDO
       WRITE(*,*) 'Ignoring unrecognized variable: ', VARN(1:NVARN)
       GO TO 10
C
 20    CONTINUE
       DO IC = 1, NCTOT
         IF(INDEX(CONNAM(IC),CONN(1:NCONN)).NE.0) GO TO 25
       ENDDO
       WRITE(*,*) 'Ignoring unrecognized constraint: ', CONN(1:NCONN)
       GO TO 10
C
 25    CONTINUE
       READ(LINE(KEQU+1:80),*,ERR=80) CONV
C
       ICON(IV,IR) = IC
       CONVAL(IC,IR) = CONV
C
      ELSEIF(KARR.EQ.0 .AND. KEQU.NE.0) THEN
C----- run case parameter data line
       PARN = LINE(1:KEQU-1)
       CALL STRIP(PARN,NPARN)
       DO IP = 1, NPTOT
         IF(INDEX(PARNAM(IP),PARN(1:NPARN)).NE.0) GO TO 30
       ENDDO
       WRITE(*,*) 'Ignoring unrecognized parameter: ', PARN(1:NPARN)
       GO TO 10
C
 30    CONTINUE
       READ(LINE(KEQU+1:80),*,ERR=80) PARV
C
       PARVAL(IP,IR) = PARV
C
      ENDIF
C
C---- keep reading lines
      GO TO 10
C
C==============================================================
C
 50   CONTINUE
      CLOSE(LU)
      ERROR = .FALSE.
      RETURN
C
 80   CONTINUE
      CALL STRIP(FNAME,NFN)
      CALL STRIP(LINE ,NLI)
      WRITE(*,8000) FNAME(1:NFN), ILINE, LINE(1:NLI)
 8000 FORMAT(/' Run case file  ',A,'  read error on line', I4,':',A)
      CLOSE(LU)
      ERROR = .TRUE.
      NRUN = 0
      RETURN
C
 90   CONTINUE
      CALL STRIP(FNAME,NFN)
      WRITE(*,9000) FNAME(1:NFN)
 9000 FORMAT(/' Run case file  ',A,'  open error')
      ERROR = .TRUE.
      RETURN
      END ! RUNGET



      SUBROUTINE RUNSAV(LU)
      INCLUDE 'AVL.INC'
C
      DO IR = 1, NRUN
        WRITE(LU,1010) IR, RTITLE(IR)
        DO IV = 1, NVTOT
          IC = ICON(IV,IR)
          WRITE(LU,1050) VARNAM(IV), CONNAM(IC), CONVAL(IC,IR)
        ENDDO
C
        WRITE(LU,*)
C
        DO IP = 1, NPTOT
          WRITE(LU,1080) PARNAM(IP), PARVAL(IP,IR), PARUNCH(IP)
        ENDDO
      ENDDO
C
 1010 FORMAT(/' ---------------------------------------------'
     &       /' Run case', I3,':  ', A /)
 1050 FORMAT(1X,A,' ->  ', A, '=', G14.6, 1X, A)
 1080 FORMAT(1X,A,'=', G14.6, 1X, A)
C
      RETURN
      END ! RUNSAV




      SUBROUTINE EXEC(NITER,INFO,IR)
C---------------------------------------------------
C     Solves for the flow condition specified by 
C     the global operating parameters:
C
C       CONVAL(ICALFA)     alpha (deg)
C       CONVAL(ICBETA)     beta  (deg)
C       CONVAL(ICROTX)     roll_rate * Bref / 2V
C       CONVAL(ICROTY)    pitch_rate * Cref / 2V
C       CONVAL(ICROTZ)      yaw_rate * Bref / 2V
C        .
C        .
C
C---------------------------------------------------
      INCLUDE 'AVL.INC'
      REAL VSYS(IVMAX,IVMAX), VRES(IVMAX), DDC(NDMAX), WORK(IVMAX)
      INTEGER IVSYS(IVMAX)
C
C---- convergence epsilon, max angle limit (radians)
      DATA EPS, DMAX / 0.00002, 1.0 /
C
      IF(LNASA_SA) THEN
C----- NASA Std. Stability axes, X fwd, Z down
       DIR = -1.0
      ELSE
C----- Geometric Stability axes, X aft, Z up
       DIR =  1.0
      ENDIF
C
      XYZREF(1) = PARVAL(IPXCG,IR)
      XYZREF(2) = PARVAL(IPYCG,IR)
      XYZREF(3) = PARVAL(IPZCG,IR)
C
      CDREF = PARVAL(IPCD0,IR)
C
      MACH = PARVAL(IPMACH,IR)
C
      IF(MACH.NE.AMACH) THEN
C----- new Mach number invalidates close to everything that's stored
       LAIC = .FALSE.
       LSRD = .FALSE.
       LSOL = .FALSE.
       LSEN = .FALSE.
      ENDIF
C
C---- set, factor AIC matrix and induced-velocity matrix (if they don't exist)
      CALL SETUP
C
      IF(NITER.GT.0) THEN
C----- might as well directly set operating variables if they are known
       IF(ICON(IVALFA,IR).EQ.ICALFA) ALFA    = CONVAL(ICALFA,IR)*DTR
       IF(ICON(IVBETA,IR).EQ.ICBETA) BETA    = CONVAL(ICBETA,IR)*DTR
       IF(ICON(IVROTX,IR).EQ.ICROTX) WROT(1) = CONVAL(ICROTX,IR)*2./BREF
       IF(ICON(IVROTY,IR).EQ.ICROTY) WROT(2) = CONVAL(ICROTY,IR)*2./CREF
       IF(ICON(IVROTZ,IR).EQ.ICROTZ) WROT(3) = CONVAL(ICROTZ,IR)*2./BREF
      ENDIF
C
C----- set GAM_U
ccc       WRITE(*,*) ' Solving for unit-freestream vortex circulations...'
      CALL GUCALC
C
C-------------------------------------------------------------
C---- calculate initial operating state
C
C---- set VINF() vector from initial ALFA,BETA
      CALL VINFAB
C
      IF(NCONTROL.GT.0) THEN
C----- set GAM_D
ccc       WRITE(*,*) ' Solving for vortex control-var sensitivities...'
       CALL GDCALC(NCONTROL,LCONDEF,ENC_D,GAM_D)
      ENDIF
C
      IF(NDESIGN.GT.0) THEN
C----- set GAM_G
ccc       WRITE(*,*) ' Solving for vortex  design-var sensitivities...'
       CALL GDCALC(NDESIGN ,LDESDEF,ENC_G,GAM_G)
      ENDIF
C
C---- sum AIC matrices to get GAM,SRC,DBL
      CALL GAMSUM
C
C---- sum AIC matrices to get WC,WV
      CALL VELSUM
C
C---- compute forces
      CALL AERO
C
C---- Newton loop for operating variables
      DO 190 ITER = 1, NITER
C
        IF(LSA_RATES) THEN
C-------- rates specified in NASA stability-axes, transform to body axes
          CA = COS(ALFA)
          SA = SIN(ALFA)
          CA_A = -SA
          SA_A =  CA
         ELSE
C-------- rates specified in body-axes, no transformation
          CA = 1.0
          SA = 0.0
          CA_A = 0.
          SA_A = 0.
        ENDIF
C
        DO K=1, IVMAX
          DO L=1, IVMAX
            VSYS(K,L) = 0.
          ENDDO
        ENDDO
C
C------ set up Newton system:  set constraints for all parameters
        DO 100 IV = 1, NVTOT
C
C-------- set index and value of constraint for this parameter
          IC = ICON(IV,IR)
C
C------------------------------------
          IF    (IC.EQ.ICALFA) THEN
           VRES(IV) = ALFA - CONVAL(IC,IR)*DTR
           VSYS(IV,IVALFA) = 1.0
C
C------------------------------------
          ELSEIF(IC.EQ.ICBETA) THEN
           VRES(IV) = BETA - CONVAL(IC,IR)*DTR
           VSYS(IV,IVBETA) = 1.0
C
C------------------------------------
          ELSEIF(IC.EQ.ICROTX) THEN
           VRES(IV) = (WROT(1)*CA + WROT(3)*SA)*DIR
     &              - CONVAL(IC,IR)*2.0/BREF
           VSYS(IV,IVROTX) = CA*DIR
           VSYS(IV,IVROTZ) = SA*DIR
           VSYS(IV,IVALFA) = (WROT(1)*CA_A + WROT(3)*SA_A)*DIR
C
C------------------------------------
          ELSEIF(IC.EQ.ICROTY) THEN
           VRES(IV) = WROT(2)
     &              - CONVAL(IC,IR)*2.0/CREF
           VSYS(IV,IVROTY) = 1.0
C
C------------------------------------
          ELSEIF(IC.EQ.ICROTZ) THEN
           VRES(IV) = (WROT(3)*CA - WROT(1)*SA)*DIR
     &              - CONVAL(IC,IR)*2.0/BREF
           VSYS(IV,IVROTX) = -SA*DIR
           VSYS(IV,IVROTZ) =  CA*DIR
           VSYS(IV,IVALFA) = (WROT(3)*CA_A - WROT(1)*SA_A)*DIR
C
C------------------------------------
          ELSEIF(IC.EQ.ICCL  ) THEN
           VRES(IV) = CLTOT - CONVAL(IC,IR)
           VSYS(IV,IVALFA) = CLTOT_U(1)*VINF_A(1)
     &                     + CLTOT_U(2)*VINF_A(2)
     &                     + CLTOT_U(3)*VINF_A(3) + CLTOT_A
           VSYS(IV,IVBETA) = CLTOT_U(1)*VINF_B(1)
     &                     + CLTOT_U(2)*VINF_B(2)
     &                     + CLTOT_U(3)*VINF_B(3)
           VSYS(IV,IVROTX) = CLTOT_U(4)
           VSYS(IV,IVROTY) = CLTOT_U(5)
           VSYS(IV,IVROTZ) = CLTOT_U(6)

C
           DO N = 1, NCONTROL
             NV = IVTOT + N
             VSYS(IV,NV) = CLTOT_D(N)
           ENDDO
C
C------------------------------------
          ELSEIF(IC.EQ.ICCY  ) THEN
           VRES(IV) = CYTOT - CONVAL(IC,IR)
           VSYS(IV,IVALFA) = CYTOT_U(1)*VINF_A(1)
     &                     + CYTOT_U(2)*VINF_A(2)
     &                     + CYTOT_U(3)*VINF_A(3)
           VSYS(IV,IVBETA) = CYTOT_U(1)*VINF_B(1)
     &                     + CYTOT_U(2)*VINF_B(2)
     &                     + CYTOT_U(3)*VINF_B(3)
           VSYS(IV,IVROTX) = CYTOT_U(4)
           VSYS(IV,IVROTY) = CYTOT_U(5)
           VSYS(IV,IVROTZ) = CYTOT_U(6)
C
           DO N = 1, NCONTROL
             NV = IVTOT + N
             VSYS(IV,NV) = CYTOT_D(N)
           ENDDO
C
C------------------------------------
          ELSEIF(IC.EQ.ICMOMX) THEN
           VRES(IV) = (CRTOT*CA + CNTOT*SA)*DIR - CONVAL(IC,IR)
           VSYS(IV,IVALFA) = ( CRTOT_U(1)*VINF_A(1)
     &                        +CRTOT_U(2)*VINF_A(2)
     &                        +CRTOT_U(3)*VINF_A(3))*CA*DIR
     &                     + ( CNTOT_U(1)*VINF_A(1)
     &                        +CNTOT_U(2)*VINF_A(2)
     &                        +CNTOT_U(3)*VINF_A(3))*SA*DIR
     &                     + (CRTOT*CA_A + CNTOT*SA_A)*DIR
           VSYS(IV,IVBETA) = ( CRTOT_U(1)*VINF_B(1)
     &                        +CRTOT_U(2)*VINF_B(2)
     &                        +CRTOT_U(3)*VINF_B(3))*CA*DIR
     &                     + ( CNTOT_U(1)*VINF_B(1)
     &                        +CNTOT_U(2)*VINF_B(2)
     &                        +CNTOT_U(3)*VINF_B(3))*SA*DIR
           VSYS(IV,IVROTX) = (CRTOT_U(4)*CA + CNTOT_U(4)*SA)*DIR
           VSYS(IV,IVROTY) = (CRTOT_U(5)*CA + CNTOT_U(5)*SA)*DIR
           VSYS(IV,IVROTZ) = (CRTOT_U(6)*CA + CNTOT_U(6)*SA)*DIR
C
           DO N = 1, NCONTROL
             NV = IVTOT + N
             VSYS(IV,NV) = (CRTOT_D(N)*CA + CNTOT_D(N)*SA)*DIR
           ENDDO
C
C------------------------------------
          ELSEIF(IC.EQ.ICMOMY) THEN
           VRES(IV) = CMTOT - CONVAL(IC,IR)
           VSYS(IV,IVALFA) = CMTOT_U(1)*VINF_A(1)
     &                     + CMTOT_U(2)*VINF_A(2)
     &                     + CMTOT_U(3)*VINF_A(3)
           VSYS(IV,IVBETA) = CMTOT_U(1)*VINF_B(1)
     &                     + CMTOT_U(2)*VINF_B(2)
     &                     + CMTOT_U(3)*VINF_B(3)
           VSYS(IV,IVROTX) = CMTOT_U(4)
           VSYS(IV,IVROTY) = CMTOT_U(5)
           VSYS(IV,IVROTZ) = CMTOT_U(6)
C
           DO N = 1, NCONTROL
             NV = IVTOT + N
             VSYS(IV,NV) = CMTOT_D(N)
           ENDDO
C
C------------------------------------
          ELSEIF(IC.EQ.ICMOMZ) THEN
           VRES(IV) = (CNTOT*CA - CRTOT*SA)*DIR - CONVAL(IC,IR)
           VSYS(IV,IVALFA) = ( CNTOT_U(1)*VINF_A(1)
     &                        +CNTOT_U(2)*VINF_A(2)
     &                        +CNTOT_U(3)*VINF_A(3))*CA*DIR
     &                     - ( CRTOT_U(1)*VINF_A(1)
     &                        +CRTOT_U(2)*VINF_A(2)
     &                        +CRTOT_U(3)*VINF_A(3))*SA*DIR
     &                     + (CNTOT*CA_A - CRTOT*SA_A)*DIR
           VSYS(IV,IVBETA) = ( CNTOT_U(1)*VINF_B(1)
     &                        +CNTOT_U(2)*VINF_B(2)
     &                        +CNTOT_U(3)*VINF_B(3))*CA*DIR
     &                     - ( CRTOT_U(1)*VINF_B(1)
     &                        +CRTOT_U(2)*VINF_B(2)
     &                        +CRTOT_U(3)*VINF_B(3))*SA*DIR
           VSYS(IV,IVROTX) = (CNTOT_U(4)*CA - CRTOT_U(4)*SA)*DIR
           VSYS(IV,IVROTY) = (CNTOT_U(5)*CA - CRTOT_U(5)*SA)*DIR
           VSYS(IV,IVROTZ) = (CNTOT_U(6)*CA - CRTOT_U(6)*SA)*DIR
C
           DO N = 1, NCONTROL
             NV = IVTOT + N
             VSYS(IV,NV) = (CNTOT_D(N)*CA - CRTOT_D(N)*SA)*DIR
           ENDDO
C
C------------------------------------
          ELSE
           DO N = 1, NCONTROL
             ICCON = ICTOT + N
             IVCON = IVTOT + N
             IF(IC.EQ.ICCON) THEN
              VRES(IV) = DELCON(N) - CONVAL(ICCON,IR)
              VSYS(IV,IVCON) = 1.0
              GO TO 100
             ENDIF
           ENDDO
C
           WRITE(*,*) '? Illegal constraint index: ', IC
          ENDIF
C
 100    CONTINUE
C

c        write(*,*)
c        do k = 1, nvtot
c          write(*,'(1x,40f9.4)') (vsys(k,l), l=1, nvtot), vres(k)
c        enddo
c        write(*,*)
C
C------ LU-factor,  and back-substitute RHS
        CALL LUDCMP(IVMAX,NVTOT,VSYS,IVSYS,WORK)
        CALL BAKSUB(IVMAX,NVTOT,VSYS,IVSYS,VRES)
C
C------ set Newton deltas
        DAL = -VRES(IVALFA)
        DBE = -VRES(IVBETA)
        DWX = -VRES(IVROTX)
        DWY = -VRES(IVROTY)
        DWZ = -VRES(IVROTZ)
        DO N = 1, NCONTROL
          IV = IVTOT + N
          DDC(N) = -VRES(IV)
        ENDDO
C
        IF(INFO .GE. 1) THEN
C------- display Newton deltas
         IF(ITER.EQ.1) THEN
          WRITE(*,*)
          WRITE(*,1902) 'iter',
     &            ' d(alpha)  ',
     &            ' d(beta)   ',
     &            ' d(pb/2V)  ',
     &            ' d(qc/2V)  ',
     &            ' d(rb/2V)  ',
     &            (DNAME(K), K=1, NCONTROL)
 1902     FORMAT(1X,A4,5A11,1X,30A11)
         ENDIF
         WRITE(*,1905) ITER, 
     &                 DAL/DTR, DBE/DTR, 
     &                 DWX*BREF/2.0, DWY*CREF/2.0, DWZ*BREF/2.0,
     &                 (DDC(K), K=1, NCONTROL)
 1905    FORMAT(1X,I3,40E11.3)
        ENDIF
C
C------ limits on angles and rates
        DMAXA = DMAX
        DMAXR = 5.0*DMAX/BREF
C
C------ if changes are too big, configuration is probably untrimmable
        IF(ABS(ALFA+DAL).GT.DMAXA) THEN
         WRITE(*,*) 'Cannot trim.  Alpha too large.  a =',(ALFA+DAL)/DTR
         RETURN
        ENDIF
C
        IF(ABS(BETA+DBE).GT.DMAXA) THEN
         WRITE(*,*) 'Cannot trim.  Beta too large.  b =',(BETA+DBE)/DTR
         RETURN
        ENDIF
C
        IF(ABS(WROT(1)+DWX).GT.DMAXR) THEN
         WRITE(*,*) 'Cannot trim.  Roll rate too large.  pb/2V =', 
     &               (WROT(1)+DWX)*BREF*0.5
         RETURN
        ENDIF
C
        IF(ABS(WROT(2)+DWY).GT.DMAXR) THEN
         WRITE(*,*) 'Cannot trim.  Pitch rate too large.  qc/2V =',
     &               (WROT(2)+DWY)*CREF*0.5
         RETURN
        ENDIF
C
        IF(ABS(WROT(3)+DWZ).GT.DMAXR) THEN
         WRITE(*,*) 'Cannot trim.  Yaw rate too large.  rb/2V =',
     &               (WROT(3)+DWZ)*BREF*0.5
         RETURN
        ENDIF
C
C------ update
        ALFA  = ALFA  + DAL
        BETA  = BETA  + DBE
        WROT(1) = WROT(1) + DWX
        WROT(2) = WROT(2) + DWY
        WROT(3) = WROT(3) + DWZ
        DO K = 1, NCONTROL
          DELCON(K) = DELCON(K) + DDC(K)
        ENDDO
C
C
C------ set VINF() vector from new ALFA,BETA
        CALL VINFAB
C
        IF(NCONTROL.GT.0) THEN
C------- set new GAM_D
         CALL GDCALC(NCONTROL,LCONDEF,ENC_D,GAM_D)
        ENDIF
C
        IF(NDESIGN.GT.0) THEN
C------- set new GAM_G
         CALL GDCALC(NDESIGN ,LDESDEF,ENC_G,GAM_G)
        ENDIF
C
C------ sum AIC matrices to get GAM,SRC,DBL
        CALL GAMSUM
C
C------ sum AIC matrices to get WC,WV
        CALL VELSUM
C
C------ compute forces
        CALL AERO
C
C
C------ convergence check
        DELMAX = MAX( ABS(DAL), 
     &                ABS(DBE),
     &                ABS(DWX*BREF/2.0),
     &                ABS(DWY*CREF/2.0),
     &                ABS(DWZ*BREF/2.0) )
        DO K = 1, NCONTROL
          DELMAX = MAX( DELMAX , ABS(DDC(K)) )
        ENDDO
C
        IF(DELMAX.LT.EPS) THEN
         LSOL = .TRUE.
C------- mark trim case as being converged
         ITRIM(IR) = IABS(ITRIM(IR))
         GO TO 191
        ENDIF
C
 190  CONTINUE
      IF(NITER.GT.0) THEN
       WRITE(*,*) 'Trim convergence failed'
       LSOL = .FALSE.
       RETURN
      ENDIF
C
 191  CONTINUE
      PARVAL(IPALFA,IR) = ALFA/DTR
      PARVAL(IPBETA,IR) = BETA/DTR
      PARVAL(IPROTX,IR) = WROT(1)*0.5*BREF
      PARVAL(IPROTY,IR) = WROT(2)*0.5*CREF
      PARVAL(IPROTZ,IR) = WROT(3)*0.5*BREF
      PARVAL(IPCL  ,IR) = CLTOT
C
      LSEN = .TRUE.
      RETURN
C
      END ! EXEC

